#include "structMesh.hpp"
#include "macros.hpp"
#include "geometry.hpp"

namespace hms {

Matrix2Xs checkVertexSizes(
	Matrix2Xs&& vertices,
	[[maybe_unused]] const Array2i& size
){
	assert( vertices.cols() == (size + 1).prod() &&
		"Number of vertices and size passed to StructMesh ctor do not match!"
	);
	return vertices;
}


StructMesh::StructMesh(
	Matrix2Xs&& vertices,
	Index nCellsX,
	Index nCellsY
) :
	StructMesh( std::move(vertices), Array2i{nCellsX, nCellsY} )
{}

StructMesh::StructMesh(
	Matrix2Xs&& vertices,
	Array2i size
) :
	m_vertices { checkVertexSizes(std::move(vertices), size) },
	m_size     { std::move(size) }
{
	m_edges.storeData();
	m_cells.storeData();
}


auto StructMesh::vertices() const -> const Matrix2Xs& {
	return m_vertices;
}

auto StructMesh::sizeImp() const -> SizeType { return m_size; }

auto StructMesh::nCellsXImp() const -> Index { return m_size.x(); }
auto StructMesh::nCellsYImp() const -> Index { return m_size.y(); }

auto StructMesh::vertexImp(Index i) const -> Traits::VertexType {
	return m_vertices.col(i);
}

auto StructMesh::domainAreaImp() const -> scalar {
	/* use shoelace formula with boundary vertices */
	Index
		nx { this->nCellsX() + 1 },
		ny { this->nCellsY() + 1 },
		nBoundaryVerts { nx + ny-2 + nx + ny-2 };

	/* we don't **have** to allocate an and copy the indices in,
	 * but it's a bit simpler */
	ArrayXi boundaryVerts (nBoundaryVerts);

	/* vertex indices */
	/* bottom boundary, left to right, with corners */
	boundaryVerts.head(nx).setLinSpaced(0, nx-1);
	/* right boundary, bottom to top, without corners */
	boundaryVerts.segment(nx, ny-2).setLinSpaced( 2*nx - 1, nx * (ny-1) - 1 );
	/* top boundary, right to left, with corners */
	boundaryVerts.segment(nx + ny-2, nx).setLinSpaced( nx * ny - 1, nx * (ny-1) );
	/* left boundary, top to bottom, with corners */
	boundaryVerts.tail(ny-2).setLinSpaced( nx * (ny-2), nx );

	return verticesToArea( *this, boundaryVerts );
	// scalar area {0};

	// Index
	// 	nxOuter { nx+2 },
	// 	domainBeg { nxOuter+1 };

	// OMP_PRAGMA(parallel for reduction (+:area))
	// for ( int j=0; j<ny; ++j )
	// 	area += cells().areas().segment( domainBeg + j*nxOuter, nx ).sum();

	// return area;
}

}